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^ Abstract 

cn 

The aim of this paper is to devise a turbulence model for the particle method 
^ Smoothed Particle Hydrodynamics (SPH) which makes few assumptions, conserves 

T— I 

linear and angular momentum, satisfies a discrete version of Kelvin's circulation 

theorem, and is computationally efficient. These aims are achieved. Furthermore, 

t> the results from the model are in good agreement with the experimental and com- 

^ putational results of Clercx and Heijst for two dimensional turbulence inside a box 

H 

^ with no-slip walls. The model is based on a Lagrangian similar to that used for the 

Lagrangian averaged Navier Stokes (LANS) turbulence model, but with a differ- 
ent smoothed velocity. The smoothed velocity preserves the shape of the spectrum 
of the unsmoothed velocity, but reduces the magnitude for short length scales by 
an amount which depends on a parameter e. We call this the SPH-e model. The 
effectiveness of the model is indicated by the fact that the second order velocity 
correlation function calculated using the smoothed velocity and a coarse resolution. 
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is in good agreement with a calculation using a resolution which is finer by a factor 
2, and therefore requires 8 times as much work to integrate to the same time. 
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1 Introduction 



The subject of this paper is a turbulence model for Smoothed Particle Hydrodynamics 
(SPH, for a review see Monaghan 2005). The SPH method was developed for astrophysical 
applications but has been progressively extended to problems involving incompressible 
fluids by using either a shghtly compressible model of the fluid (Monaghan 1994), or by 
algorithms designed to solve the full incompressible equations (Cummins and Rudman 
1999, Hu and Adams 2007). In this paper we use the slightly compressible model. 

Many of the incompressible flow problems to which SPH has been applied involve 
the generation of turbulence. A common example is a dam break where the flow of the 
escaping fluid is laminar until it hits a downstream wall and forms a breaking wave. 
Another example is the sloshing of a fluid in an oscillating tank where breaking waves are 
created each cycle of sloshing. Standard turbulence models such as the k — e model have 
been apphed to the SPH simulation of dam break problems (Violeau and Issa 2007) with 
modest success. In this paper we make use of the ideas associated with the Lagrangian 
averaged Navier Stokes equations (LANS) explored by Holm and his colleagues (Holm 
1999, 2002, Chen et al. 1999, Cheskidov et al. 2005, Guerts and Holm 2006), Mohscni et 
al. 2003, and for extensive references see Graham et al. (2007) and Lunasin et al. (2007, 
2008). The basic idea is to determine a smoothed velocity v by a hnear operation on 
the un-smoothed velocity v, and then determine the Eulerian equations of motion from 
Lagrange's equations for the SPH particles using a Lagrangian where the kinetic energy 
per unit mass is |v • v. 

The average motion of the fluid is determined by v and in the SPH formulation the 
particles are moved with this velocity. The resulting acceleration equation contains extra 
terms which represent the effective stresses induced by the smoothing. Once the form of 
the smoothing is chosen these stresses are determined. The equations of motion conserve 
energy, linear and angular momentum (in the absence of rigid boundaries and external 
forces), and they satisfy a discrete version of the circulation theorem. 

To complete the model we add a viscous dissipation term. We use one of the 
standard SPH viscosity terms which has been tested for Couette flow (Monaghan 2006), 
and for spin down in a rotating cyhnder (Monaghan 2005, Monaghan and Kajtar 2009). 
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We simulate he boundary conditions associated with fixed or moving bodies by placing 
boundary particles on the boundaries. These boundary particles exert forces on the fiuid 
SPH particles. This technique, which is similar to the Immersed Boundary Method of 
Peskin (1977, 2002), has been shown to give good results for a wide range of problems 
(Monaghan and Kajtar 2009). The effects of the turbulence on processes such as thermal 
transport can be estimated using the differences between the smoothed and unsmoothed 
velocities to determine diffusion coefficients. However, we do not study those processes in 
this paper. 

We apply the SPH-e equations to turbulence in a two dimensional box with no-slip 
boundary conditions. This system has been studied both experimentally and numerically 
by Clerx et al. (1999, 2000), Massen et al. 2002) and is particularly interesting because 
it shows that the dynamics with no-slip or stress-free boundary conditions (relevant to 
many turbulence problems in nature), differs significantly from that with periodic condi- 
tions. This includes the formation of vortices, the spectrum, and changes in the angular 
momentum. The two dimensional problem also has the advantage that the computational 
demands are very much less than for a three dimensional simulation. 

The plan of the paper is to discuss the smoothing, and show how it can be formulated 
so that the linear and angular momentum are conserved in the absence of boundaries and 
external body forces. We then derive the Euler equations of motion from a Lagrangian 
and complete the SPH-e equations by adding the boundary force and viscous terms. We 
discuss the conditions for the kinetic energy to be positive definite and describe the time 
stepping scheme. We apply the SPH-e equations to predict the properties of decaying 
turbulence initialised by a set of gaussian vortices and compare them with those of Clercx 
et. al. (1999, 2000), and Massen et al. (2002). 

2 Smoothing 

The typical LANS model uses a smoothed velocity v defined in terms of the the disordered 
velocity v by 




(2.1) 
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where the integration is over the region occupied by the fluid, and G is a kernel which 
satisfies 

^G(|r'-r|,£)dr' = 1, (2.2) 



and is a member of a sequence of functions which tends to the S function in the limit 
where £ — > 0. A typical example is a Gaussian, though in practice we use smooth functions 
that have compact support. The length scale i determines the characteristic width of the 
kernel. For the reader famihar with SPH it is useful to note that the kernel G has the 
same properties as the kernel W used in SPH, and i has the same significance as the 
length scale h used with the SPH kernel W. In the applications to be described here we 
replace G by the W used in the SPH fluid dynamics and we replace £ by h. 

It is common practice (e.g. Chen et al. (1999), Lunasin et al. 2007) to replace the 
integral smoothing by the implicit differential smoothing 

v{r) = {l-a'V')9{r), (2.3) 

where is the Laplacian operator, and a has the dimensions of length and determines 
the length scale below which velocity variations arc smoothed. If v and v arc expanded 
in a Fourier series with coefficients of the Fourier term e^^'^, Ck and Ck respectively, then 
equation (2.1) gives 

Ck = G{k)Ck (2.4) 

where G{k) is the Fourier transform of G. If G is a smooth enough function the Fourier 
transform decreases rapidly with increasing k and Ck <^ Ck The differential smoothing 
(2.3) behaves similarly and leads to 

A disadvantage with (2.3), and with the smoothing discussed by Monaghan (2002) for 
a preliminary form of SPH turbulence, is that it is implicit and can only be solved by 
iteration when the domain is complicated and Fourier spectral methods cannot be used. 
In three dimensions the large number of iterations required may then make the implicit 
smoothing impracticable. The implicit SPH smoothing (Monaghan 2002) was too slow 
even for two dimensional problems because of the large number of iterations required. 
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and the implicit smoothing (2.3) greatly increases the computational time of LANS-a 
ocean models (Hecht et al. 2008). We therefore consider a different smoothing which 
can be converted to the SPH formulation easily, and involves neghgible extra work. The 
smoothing is defined by 

v(r) = v(r) + e J (v(r') - v(r))G'(|r' - r|, £)dr', (2.6) 

where e is a constant and < e < 1. In this case the Fourier coefficients satisfy 

dk^Ck{l + e{G{k)-l)), (2.7) 

and Cfc — >• (1 — e)Cfc as /c — oo. Provided (1 — e) ^ 1 and < e < 1 then -C Cfe as A; ^ 
oo. The smoothed velocity spectrum therefore retains the same form as the unsmoothed 
spectrum, though it is reduced in magnitude. Of course, in a dynamical calculation, the 
spectrum of a quantity such as the kinetic energy depends on the smoothing in a more 
comphcated way. The standard LANS-a model changes the spectrum for high k as in 
(2.5). In the calculations we describe in this paper, the typical value of e is 0.8, so that 
the high order Fourier coefficients are reduced by a factor 0.2 every step. We refer to this 
model as the SPH-e turbulence model. 

The SPH equivalent of (2.6) is the XSPH equation (Monaghan 1989) which can 
be obtained by using the SPH procedure for converting integrals into summations over 
particles. This procedure begins with the integral interpolant which estimates scalar, 
vector or tensor function T at any point by T/ defined by 

r,(r) = jT{r')W{\r'-r\,h)dr', (2.8) 

where W is a smoothing kernel normalized as for G in (2.2), and /i is a length scale 
associated with the kernel. We use kernels with compact support which vanish for points 
separated by more than 2h. 

This integral can be approximated by summing over particles (the errors involved 
with this are discussed by Monaghan 2005) to give 



r(r) = J]^W((|r-r,U), (2.9) 
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where m^, Tj, and pb are the mass, value of T and the density of particle b. The subscript 
I in (2.8) has been dropped for simplicity. An example of this interpolation is the density 
which can be estimated anywhere by 

p(r) = J]m6iy((|r-r6U). (2.10) 

b 

For other examples see Monaghan(2005). The summation is over all particles but the only 
contributions to the value of T at particle a are those within 2h. Efficient methods of 
finding these particles means that the work done to find T for all particles is proportional 
to the number of particles. Spatial derivatives of quantities such as T(r) can be obtained 
at the position of any particle by differentiating (2.9) analytically and evaluating the 
resulting expression at the position of that particle. 

Applying this procedure to (2.6) the smoothed velocity for particle a becomes 

^^a = v„ + e V— (vfe-Va)^^^, (2.11) 
b 

Where Gab denotes G(|ra — r^l, £). However, because the centre of mass should move with 
constant velocity when the fluid is isolated and free of external forces, we require that 

^ mar a = XI = XI "^oVa = ^ ^ ^' ^'^'^'^^ 

a a a a 

where V is the constant momentum.The second equality follows from the fact that ele- 
ments of fluid move with the smoothed velocity. The last equality follows from the fact 
that we use a Lagrangian. However, to satisfy the third equality the summation term in 
(2.11) must be anti-symmetric in a and h when multiplied by ma- To this end we replace 
Pb by a symmetric density pab and £ by a symmetric form for example — \{^a + h)-i 
and write 

Gab ^ K{\ra-rb\,£ab) ^2 13) 

Pab ^ibPab 

where d is the number of dimensions. In compressible gas dynamics £^ oc 1/pa which 
means that £ is proportional to the local average particle spacing. The same principle 
applies also to the nearly incompressible fluid where the changes in density and hence i 
are small. We can now define the symmetric density by the condition 

KbPab = M, (2.14) 



where M is a mass which is held constant (suggested by Daniel Price). This mass is 
determined by the initial state where p and I are the same for each particle. The final 
form of the smoothing is then 

= Va + e ^(Vfe - Y^)Kab, (2.15) 
6 

In the absence of external forces and fixed boundaries the angular momentum is 
also conserved because X^a^^a^o x Va = 0, and then 

^ ^ m„r„ X v„ = ^ ma X v„ + r„ X ^ j = ^ mar„ x ^ = 0, (2.16) 

a a ^ 'a 

where the last equality follows from the invariance properties of the SPH-e equations (see 
below) . 

In addition to ensuring conservation of linear and angular momentum this choice of 
'pah means that wc do not have to consider the change of the density term in (2.15) when 
working out the Lagrange equations. A further simplification, suitable for the nearly 
incompressible fiuid, is to assume tab is also constant when working out the Lagrange's 
equations. 



3 Lagrange's Equations 

The Euler equations for the fiuid are the Lagrange equations obtained from the SPH 

equivalent of Eckart's Lagrangian (Eckart 1960). 

= X^J^ifo • Vft - M(pb,Sb)^ , (3.1) 

where ^(pb, s) is the internal energy per unit mass associated with the equation of state. 
This energy depends, in general, on the density p and the entropy per unit mass s which, 
for the non dissipative Lagrangian formulation, does not vary with time. In this paper 
we assume each particle has the same entropy. The invariant energy (in the absence of 
viscous dissipation and external or boundary forces) is 

£; = ^7n6 Qvb-V6 + M(p6,S6)^ . (3.2) 
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Lagranges equations for particle c are 

d f dL\ dL 



dt \ d\r / dr, 



(3.3) 



The Lagrangian must therefore be written in terms of the v and r before the Lagrange 
equations are worked out. In order to do this it is convenient to write (2.15) in the form 

Pa = rUaVa = ^afoVfo, (3.4) 
b 

where 

Dab = 5ab{ma - ^TajKaj) + e(l - 5ab)(7abKab, (3.5) 

Cab = TnaTrih/M and 5ab is the Kronecker delta. Note that Dab is symmetric and the matrix 
D with elements Dab is a square, symmetric matrix. We can then write 

b 

where {D)~^ is the ab component of the matrix inverse to D. The kinetic energy 

-Eif = ^ ^ m„Va • Vtt, (3.7) 

a 

can then be written in the form 

EK^lj2P'^-J2(^^abPb. (3.8) 



a 



Therefore 



^^^{maSacPb + 'mb6bcPa)Dj. (3.9) 
a b 

The first and second terms in the summation are identical, and equal to IrricVc, so that 
the canonical momentum is rric^c- 

The right hand side of Lagrange's equations (3.3) has two contributions, one from 
kinetic energy term and one from the elastic energy. The second leads to the usual SPH 
pressure equations (see for example Monaghan 1992, 2005). The first can be written 

BEk d 



(pD-ip) = p (^^D-^) p. (3.10) 



where p = (pi, P2, Ps, ■ ■ ■), and the subscript denotes a particle label. The row or column 
form is taken as required for the matrix products. The derivative d/dvc is to be interpreted 
in terms of a separate matrix equation for each component of Yc- Because D is a square 
symmetric matrix so is its inverse, and we can use the identity 



If this result is substituted into (3.10) and note is taken of (3.4), we can write (3.10) as 



where v = (vi, V2, V3, • • •) and in (3.12) the column vector form is used. We can then 
write 



The expression (3.5) for Dah contains two summation terms and in each case they involve 
derivatives of the kernel. However, the discussion after (2.12) shows that we need only 
consider the explicit dependence on the coordinates, and neglect the variation of lab and 
Pa}}. If the fluid is compressible then the variation of I ah involves extra spatial derivatives 
but in this paper we assume the variations in lab are negligible. The first contribution to 
Dab is 




(3.11) 




(3.12) 




(3.13) 




(3.14) 



with derivative 




(3.15) 



The contribution to the derivative of the kinetic energy is 




(3.16) 



If this expression is simplified, noting that dKab/dva 



dKab/ drf,, it becomes 




(3.17) 
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The second contribution to Dab is 

D^^ = e(l - Sab)aabKab. (3.18) 

with derivative 

^ = 6(1 - 6ab)aab^{6ac ' ^bc), (3-19) 

and the contribution to the derivative of the kinetic energy is 

^ = • ^b(l - Sa,)aa,^{Sa. - 5,.). (3.20) 

If this expression is simphfied it can be written 

dE^i^ ^ dK, 



OTc dr. 

Combining (3.17) and (3.21) gives 

OEk 1 o dK, 



-e^(T„cV„- Vc^-^. (3.21) 



J^E->«^^' (3-22) 



where Vab = Va - v^. 

The contribution to the equations of motion from the internal energy term can be 
calculated easily. The details can be found in Monaghan (2005), but the essentials are 
as follows. If there is no dissipation the first law of thermodynamics for unit mass is 
du = dp/ p^. For the slightly compressible fluid we assume 



where po is a reference density, Cs is the speed of sound, and 7 = 7 in this paper. For 
the weakly compressible case we choose Cg — lOV^ai where V^ax is the maximum fluid 
velocity, so that the density fluctuation 5p/ p ~ {vmax/csY — 0.01. If the density is 
estimated by the SPH summation 

b 

the contribution to the right hand side of Lagrange's equation for particle c 

Edui, dpb 
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becomes 



(3.26) 



Combining the previous results, Lagrange's equation for particle c are 




(3.27) 



This equation is the SPH-e Euler equation. It is Galilean invariant and invariant to 
rotations of the coordinate system. For these reasons the linear and angular momentum 
are conserved. Because the Lagrangian is not an explicit function of time the energy (3.2) 
is also conserved. The equations also satisfy a discrete version of Kelvin's circulation theo- 
rem (this follows from the same argument used by Monaghan (2005)). The term involving 
e in (3.27) is related to the velocity derivative terms in the LANS-alpha equations. 

In this paper we choose the smoothing function G to be the same as W and £ to be 
h. In that case we can replace (3.27) by 



where h is allowed to vary according to ha oc l/pj This variation is < 1% when the 
speed of sound is calculated as described after (3.23). 

In addition to the acceleration equation just derived we use a density continuity 
equation for particle c which may be deduced directly from the continuum continuity 
equation, or by differentiating (3.24) with respect to time. It is 



where v^c = — Vc, and 




(3.30) 



The acceleration equation (3.27) can also be derived using the SPH continuity equation 
as a constraint equation when deducing Lagrange's equations from the Least Action vari- 
ational principle. 




(3.28) 




(3.29) 



b 
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3.1 Boundaries and boundary forces 

In this paper the boundaries are defined by boundary force particles which exert forces 
on the fluid. This technique is similar to the immersed boundary method of Peskin (1957, 
2002) and is discussed in detail by Monaghan and Kajtar (2009). The basic idea is that 
if the interaction between the fluid and boundary particles is sufficiently smooth, and the 
forces between the particles is along their fine of centres, the force on a fluid particle is 
normal to the boundary high accuracy. Monaghan and Kajtar (2009) showed that this 
is true for both curved and planar surfaces, provided the boundary force particles were 
spaced less than half the spacing of the fiuid particles. 

When the boundary force particles are included (3.28) becomes 

t = - E-"' (§ + 5 - #) Vcl^-c. + t E f.>. (3-31) 

b ^Pbc/ k=ijes, 

where the sum over k is over the Nf, bodies, and j e Sk denotes the boundary particles 
of body k. The details of the boundary force are given in the appendix. Because the 
boundary force between a pair of particles is radial it is possible to write it in terms of a 
pair potential which can then be included in the Lagrangian. This boundary force does 
not change the conservation properties of the Lagrangian. 

The present formulation is general and allows for an arbitrary number of boundaries 
which could include immersed bodies. These bodies may move either by specifying their 
motion or by determining it from the forces exerted on the boundary force particles by the 
fluid particles. In this paper the only boundary is the fixed square boundary containing 
the fluid. 

3.2 Viscous forces 

To complete our model we need to add a viscous term. There are several SPH approxima- 
tions to the Navier-Stokes viscosity (see Monaghan 2005 for two examples) . The viscous 
equations take the form 

^ = - ^ + ^ + - VeW^,, + E E (fci - n,cVeW^,e) ■ (3.32) 
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where 

^^^^_av^^^n^c:nc_ (3.33) 

Pbc^bc 

In this expression a is a constant, Vgig is a signal velocity which, in this paper, is taken 
as the speed of sound in the equation of state, and p^c — {pb + Pc)/'^- If the density is 
constant, and the kernel 1^ is a Wendland fourth order kernel (see §3.3), it is possible to 
show that the kinematic viscosity v is then given by 

u — -acsh. (3.34) 
8 

This result is obtained by converting the SPH summations to integrals (Monaghan 2005). 
3.3 Constraints on e from the kinetic energy 

The kinetic energy should be positive definite. From (3.4) and (3.7) the kinetic energy is 

^x = ^ J]J]v,Z},,Vfc, (3.35) 



2 

a b 



which can also be written as 



ek-Iy: - E E (3-36) 

a a b 

which shows that is less than the kinetic energy using the unsmoothed velocity. This 
is to be expected because the smoothing reduces the velocity associated with short length 
scales. However, it is important to ensure that Ek is positive definite otherwise the 
energy might be reduced by increasing the velocity. Clearly this condition depends on the 
magnitude of e. In the calculations described in this paper where < e < 0.9 we always 
find Ek > 0. 

A sufficient condition on e can be found by determining the sufficient condition for 
the eigenvalues of the matrix D to be positive. This can be done using Gerschgorin's 
theorem. This simple theorem is obtained as follows. Let Dq = Aq. We normalize the 
eigenvector q by making the largest component Qk — 1- Then, from row k of the matrix 
D, we deduce 

E + ^'<^k = A, (3.37) 
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and 

Prom this inequality we infer that A > provided 

D^k>Y,\^kj\. (3.39) 

3+k 

or 

To estimate the values of e for which this inequality is true we consider equal mass SPH 
particles placed on the vertices of a grid of squares where the side of each square has 
length A. In the calculations to be described here the smoothing function is the fourth 
order Wendland kernel for two dimensions (the order here refers to the way it approaches 
zero) which has the form 

^^^'^^ = 64^^^"^/^^'^^ + ^^/^^' ^^-^^^ 

\l z < 21 and zero otherwise. The smoothing is therefore over a length 21. If £ = A 
the sufficient condition is that e < 1.04. Hi— 1.5 A the conditions is e < 0.667. If i is 
sufficiently large relative to A we can replace the summations by an integration over the 
area after subtracting off the term with j — k m. (3.6). We then find that the sufficient 
condition is 1 > 2e(l — 7A^/ (47r£^)) or, when £/A ^ oo, the sufficient condition becomes 
e < 0.5. However, as mentioned earlier, these sufficient conditions appear to be excessively 
pessimistic because a wide variety of calculations with e = 0.9 have positive Ek- 

3.4 Time Stepping 

The SPH turbulence equations were integrated using a time stepping scheme that is second 
order and based on the Verlet symplectic method though it is more complicated because 
the force depends on the velocity. The ideal form is reversible in the absence of viscosity. 
For convenience we write the equation in the form 

^ = F(r,p,v)e, (3.42) 
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f = B(r.v), (3.43) 

Vc = g(vc,rc), (3.45) 

We use the notation that A° denotes a quantity A at the beginning of the current step, 
A^^^ at the midpoint of the step, and A^ at the end of the step. The time stepping 
equations, where 5t is the time step, can then be written 

= r'^ + W, (3.46) 



r 

-c ■ 2 



P'J' = + (3.47) 

The time step can then be completed by first calculating according to 

v,^ = v° + 5^Fe(rl/^ (v^ + v°)/2), (3.48) 

after which can be calculated according to 

= g(v\r^), (3.49) 
vl = vlf' + ^5tvl. (3.50) 

This algorithm is reversible in the absence of dissipation. However, iteration is required to 
solve (3.48). We therefore compromise by replacing (v^ + v''c)/2 with vl^^ and estimate 
v^^ by 

vV^ = v° + IstF:^^ (3.51) 

The two equations (3.49) and (3.50) also require iteration. Clearly the number of iterations 
depends on the time step and e and, if either is too large, the iterations will not converge. 
The condition on the time step to guarantee the iteration will converge is easily worked 
out. However, for all the calculations in this paper where < e < 0.9, the iterations 
converge in 3 steps with an error of 10~^ with just the CFL condition St < O.bh/cg. Note 
that the function g(v, r) also contains £ or, in the present case, h. When iterating (3.49) 
and (3.50) this is replaced by the midpoint value calculated using the midpoint density. 



16 



4 Applications to turbulence in a square 2D rigid box 



We consider a two dimensional fluid inside a rigid square with side length Im and no-slip 
boundary conditions. This problem provides a good test of turbulence models because 
it has been studied in great detail by laboratory experiment and by highly accurate 
numerical simulations (see for example Clercx et al. 1999, Clercx and Heijst 2000, and 
Maassen et al. 2002). Furthermore, the computer time is greatly reduced relative to three 
dimensions and the calculations described here were run on a MacBook Pro. The results 
of Clercx et al. which we use for comparison have Reynolds numbers K in the range 1000 
to 5000. In this range the results are similar. 




Figure 1: The v field for e = 0.9 shown after 5s (~1 typical time scale). The small length scale 
structures have been replaced by larger scales. 

For the SPH simulation we use the 4th order Wcndland kernel W{r, h) identical to G 
defined by (3.41). The initial particle density was l.Olpu, where pw is the density of water, 
and the reference density for the equation of state. This produces a small background 
pressure, h was set to 1.5dp where dp is the particle spacing. The smoothing of the 
velocity has i — h so that the smoothing is over a circle of radius ~ 3dp around any 
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specified particle. 

The Reynolds number 3?, with velocity 0.12 m/s (which is the typical maximum 
velocity in the box after the few seconds), and length scale equal to half the side of 
the box, was 1000. This 3? is similar to those used by Clercx, Maasscn and van Heijst 
(1999) who considered 3? in the range 500 < 3? < 2000. With the 3? specified, and the 
kinematic viscosity given by (3.34), the value of a for a specified resolution and therefore 
specified /i, can be calculated. The resolution length for a well resolved calculation can 
be estimated in the same way as did Clercx, Maassen and van Heijst (1999) from the 
enstrophy dissipation rate C using the combination {C^jv^Y^^ « (boxlength)/A'" where 
N is the number of Chebyshev modes along a coordinate axis. They conclude that for 
3? = 1000 that should be ~ 180, or the total number of modes is 180^. We find that 
convergence is achieved for an SPH calculation with 150^ particles. 

The SPH particles were placed initially on a grid of squares then damped to equi- 
librium, after which they were set in motion with velocities specified by a 4 x 4 set of 
vortices. A similar procedure was used by Maassen et al. 2002 though they had higher 
resolution and used a 10 x 10 set of Gaussian vortices. The vortices were equi-spaced, 
with spacing 0.2, then given a random shift in both x and y of 0.02 77 where 77 is a uni- 
formly distributed random number with — 1 < 77 < 1. The sign of the rotation and hence 
vorticity due to each vortex was ± in a chess board pattern. The velocity at r due to a 
vortex at R is 



and d = 0.02. 

The Fourier transform of the vorticity produced by this velocity field decreases ex- 
ponentially according to e"^'''^/^^^. The no-slip condition requires that the fiuid velocity 
vanishes on the boundary. This was achieved by smoothing the velocity near the bound- 
aries in a similar manner to Clercx et al. (2000). The mean square speed was initially 
0.15 giving a typical time scale of0.5/.15~3. 



v(r) = fle^ X (r - R), 



(4.1) 



where is a unit vector normal to the fiuid and Q is given by 




(4.2) 
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0.2 0.4 0.6 0.8 

Figure 2: The v field for e = 0.9 shown after 15s (3 typical time scales). The small length scale 
structures have been replaced by larger scales. 

In Fig.l we show the smoothed velocity field of a set of fiuid particles with e = 0.9 
and initial spacing 1/75 at a time 5s (or ~ 1 time scales) after damping. Figure 2 shows 
the smoothed velocity field after 15s (or 5 time scales. Already much of the complex short 
length scale velocity field has transformed into larger scales which is similar to that found 
in experiments and in other simulations (Maassen et al. 2002, their figure 3). A similar 
result was found if the simulation is run without smoothing indicating that it is not due 
to the smoothing but to the normal processes that occur in two dimensional turbulence. 

4.1 Decay of kinetic energy and enstrophy with time 

In Fig. 3 we show the decay of kinetic energy (defined as in (3.7)) , for three simulations 
with initial particle spacing 1/75, 1/125 and 1/150 ande = 0. The decay is similar to the 
experimental results found by Maassen et al. (2002) for weak stratification (their figure 

9(a)). The convergence of the numerical solution is indicated by the closeness of the 
results for the particle spacing 1/125 and 1/150. We estimate the error in a simulation 
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Figure 3: The decay of kinetic energy with time for a simulation with e = is shown. The filled 
circles are for resolution 1/75, the larger open circle are for resolution 1/125 and the continuous 
curve for resolution 1/150. The small open circles mark a line with time variation 

with particle spacing dp to be oc {dpY'^. Similar results were found for the enstrophy in 
the same simulation shown in figure 4 where the hne of small open circles has the time 
variation l/t^-"^ in good agreement with the results of Maassen et al. (2002) (their figure 
9(b)). In figure 5 the time variation of the ratio of enstrophy to kinetic energy is shown. 
This is also in good agreement with Maassen et al. (2002) (see their figure 17). Figure 6 
shows the decay of the enstrophy with time when e = 0.9. 

4.2 The second order velocity correlation function 

The second order, longitudinal, velocity correlation function C2{R) is defined by 
r(J,^ //[(v(r)-v(rO)-q/gf%-i?)drdr- 

where q = r — r', 5{x) denotes a one dimensional Delta function, and dr denotes a volume 
(area in 2D) element. Depending on the case we consider the velocity may be the smoothed 
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Figure 4: The decay of Enstrophy with time with e = 0. The filled circles are for resolution 
1/75, the larger open circle are for resolution 1/125 and the continuous curve for resolution 
1/150. The small open circles mark a line with time variation 
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Figure 5: The decay of Enstrophy/ Ek with time with e = 0. The filled circles are for resolution 
1 /75, the larger open circle are for resolution 1/125 and the continuous curve for resolution 1/150. 
The small open circles mark a line with time variation . 
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Figure 6: The decay of the entrosphy for initial spacing 1/75 (filled circles), and 1/125 (large 
open circles), with e = 0.9. The line of small open circles marks 1/t^'^ 

or unsmoothed velocity. We evaluate the previous expression using SPH summations. The 
most efficient way of evaluating these summations is by binning. We consider every pair 
of particles. For particles a and h we calculate an integer k = Int{rab/^) where Int is 
the Fortran integer function and A is a suitable small fraction of the maximum value of 
R. The value of ((v^ — V;,) • rah/rabf is added to and 1 is added to N^. The value 
of C2{Rk) is then the final value of F^/N^ where Rk — kA. Because the system is not 
homogeneous or isotropic the correlation function was calculated for points within the 
square 0.3 < x < 0.7, and 0.3 < y < 0.7 so that the influence of the boundary was small. 

An example of the results found using this method are shown in Fig. 7 for the case 
of no smoothing and initial particle spacings of 1/75, 1/125 and 1/150 at a time of 20s 
after the velocities were initialised. The values of C2 were scaled by the average of C2. 
All three correlation functions agree closely until a radius of approximately 3/75 which 
is close to 2h for the lowest resolution. At this length scale the correlation function with 
resolution 1/75 is larger, indicating more disorder than for the higher resolutions. The 
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Figure 7: The second order velocity correlation functions for three simulations with e = 0, 
approximately 20s (four typical time scales) after the velocity was initialized. The dashed line 
shows the results with initial spacing 1/75, the dot-dashed results with initial spacing 1/125 and 
the filled star symbols the results with initial spacing 1/150. The straight dotted line marks the 
line B}- '^. 
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good agreement between the results for initial particle separation 1/125 and 1/150 over 
the whole domain confirms the results for the decay of kinetic energy and enstrophy. The 
radial variation of the correlation function is given approximately by R^'^ with a distortion 

at large R caused by the finite domain sampled, and possibly to the no-slip boundary 
conditions. The Batchelor theory (Batchelor 1969) of the direct cascade of Enstrophy in 
two dimensional turbulence predicts that this second order correlation function should 
vary with R as R^. 




Figure 8: The second order velocity correlation function for simulations with initial particle 
spacing 1/75 and 1/150 at a times 20s after the velocities have been initialized. The dashed 
line shows the second order correlation function for the v, and the dash-dot line shows the 
correlation function for v in a simulation where e = 0.75. Note that the correlation function 
with the smoothed velocity, and particle spacing 1/75 , is very close to that with no smoothing 
and spacing 1/150. 

In Figure 8 we show the second order velocity correlation function calculated for 
both the smoothed and unsmoothed velocities when the initial particle spacing is 1/75 

25 



and e = 0.75, and a higher resolution calculation with spacing 1/150 and no smoothing. 
The star symbols show the results when the initial particle spacing is 1/150 . These results 
show that with e = 0.75 the results for the lower resolution are in very close agreement 
with those for the higher resolution, indicating that the smoothed coarse simulation gives 
a similar velocity field to the higher resolution calculation. It is interesting to note that 
the coarse resolution with no smoothing gives a velocity correlation function between that 
for V and v calculated using smoothing. 
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Figure 9: The Chebychev spectrum of the kinetic energy at a time 40s after damping with 
initial particle spacing 1/75. The continuous line is for e = and the dashed line is for e = 0.9. 
Note that the larger e results in a significant drop in the spectrum after mode number ^ 8. The 
small circles mark the line The zig zag shape shows that the odd modes have more 

energy than the even modes. 

4.3 Chebyshev spectrum 

Because the motion in a rigid no-slip box is neither homogeneous nor isotropic Clercx and 
Heijst (2000) calculated one dimensional spectra both near the boundaries and in the cen- 
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tral regions. The spectrum shown in figure 9 was calculated following the procedure used 
by Clercx et al 1999, where the kinetic energy is expanded in Chebychev polynomials. The 
version of their procedure adopted here is the following. We use Chebyshcv polynomials 

T*{x) related to the standard Chebyshev polynomials T„ by T*{x) = T„(2a; — 1) where 
< a; < 1. If Ek (defined by (3.7) ) is expanded in these polynomials it can be written 

EK{x,y) = YmC,,T:{x)T;{y). (4.4) 

i 3 

The coefficients Cij were calculated by integrations which were evaluated using SPH 
summations. Comparison with the expansion of test functions showed that to estimate 
the coefficient C*^ to within 10% required both i and j to be < (boxwidth) / (5(ip) where 
dp is the initial particle spacing. 
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Figure 10: The Chebychev spectrum of the kinetic energy with e = at a time 40s after 
damping with initial particle spacing 1/125 (shown by the dashed fine) and 1/150 (shown by 
the continuous line). The small circles mark the line 1/n^"^. 

If this expansion is evaluated for the lines x — \ and y — \ and the results combined 
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we can define a one dimensional spectrum 

ek^Y.^:t:{x), (4.5) 

n 

where 

C:^lj2(^m + Qr.mi/2), (4.6) 

i 

where the summation has only even values of i. Some details of this calculation are given 
in the Appendix. 
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Figure 11: The Chebychev spectrum of the kinetic energy for initial spacing 1/75 (with e = 0.9) 
shown by the dashed line), and 1/150 (with e = 0.0), shown by the continuous line). The small 
circles mark the line 1/n^'^. 

Figure 9 shows the Chebyshcv spectrum when the resolution is 1/75 and e has the 
values (continuous line), and 0.9 (dashed line). As in all the Chebyshev ID spectra 
(including those of Clercx and Heijst 2000) the spectrum has a pronounced zig-zag ap- 
pearance because the modes with odd n have coefficients that are larger than those for 
even n. The smoothing reduces the amplitude of the zig-zag. The spectrum is in satisfac- 
tory agreement with that of Clercx and van Heijst (2000) (see their figures 2(a)), though 
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they calculate the spectra with higher resolution and a corresponding larger number of 
modes, and they average over two eddy turnover times. Their C* fall off typically as 1/vP 
where 2 < p < 3 , and in figure 10 the line marked by small circles shows the case p = 2.5. 
The spectrum with e = 0.9 is significantly below that for e = after mode number ~ 10. 
This mode number is equivalent to a length scale of 0.1 or 5/i. 

Figure 10 shows the Chebyshev spectrum when the initial particle spacings are 1/150 
and 1/125 and there is no smoothing. This figure shows that the spectra are very similar 
for the two resolutions although the results for low mode numbers are more jagged for the 
1/150 than for 1/125. The agreement between these results confirms the inference from 
figures 5 and 7 that, at these resolutions, convergence is achieved. 

Figure 11 shows the spectrum with initial particle spacing 1/75 and e = 0.9 and the 
spectrum with initial particle spacing 1/150 and e = 0. Although the spectrum for the 
spacing 1/150 runs to higher mode numbers than for 1/75, there is a strong suggestion 
that both spectra have a similar trend for mode numbers > 8 but the averaged smoothed 
spectrum is lower, which is in agreement with the argument given in §2. No stronger 
conclusion can be drawn because the spectra are jagged. It must be kept in mind that 
the one dimensional spectra depend on the where they are constructed. For example we 
could use lines parallel to the y axis for x equal to | and | with similar lines parallel 
to the X axis. In this case the lines are closer to the no-shp boundary and the spectrum 
changes. Clercx and Heijst (2000) consider the lines we use here as well as lines close to 
the boundaries where the spectrum falls off as 1/n^/^. 

5 Conclusions 

The results of this paper show that SPH simulations using a standard algorithm reproduce 
the results obtained by other researchers for decaying turbulence in a square box with 
no-slip boundary conditions provided the resolution length is < 1/125, in agreement with 
estimates made using the Reynolds number. These results include the time variation of 
the decay of kinetic energy, enstrophy and the form of the spectra. The turbulence model 
proposed in this paper reproduces these results but, in addition, shows that quantities 
such as the second order velocity correlation function calculated using a resolution of 
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1/150, can be reproduced with an initial particle spacing 1/75 which requires a factor 8 
less computation to integrate to a specified time. 

The new algorithm is easy to implement and requires typically only 20% more com- 
puting time than a standard integration. An attractive feature of the algorithm is that, 
like the LANS-a model, it conserves linear and angular momentum (in the absence of fixed 
bodies or boundaries and external forces), and satisfies a discrete form of Kelvin's circu- 
lation theorem. Detailed exploration of the the dynamics for other boundary conditions, 
and other boundaries, together with driven turbulence is in progress. 
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8 Appendix 



8.1 Boundary forces 

In this paper we use the following form for the force per unit mass on fluid particle c by 
boundary particle j, icj, by 

i\rc,-d\) r,' ^'-'^ 
where B{r) is the one dimensional cubic Wendland function (Wendland )which has been 
scaled so that it is only non zero in the domain < r/h < 2 where it has the form 

The parameter d is the spacing of the boundary force particles. Provided d is less than 
0.5dp, where dp is the typical fluid particle spacing, the tangential component of the 
force is < 10~^ of the normal component. Further details can be found in the paper by 
Monaghan and Kajtar (2009). 

8.2 Chebyshev spectra 

The continuum version of the SPH kinetic energy is 

j j E{x,y)dxdy = - j j p'v{r) ■ \{r)dxdy (8.3) 
JO Jo 2 Jo Jo 

The expansion of the kinetic energy per unit area E{x, y) in Chebyshev polynomials T*{q), 
such that T*{q) = T„(2g — 1), with < g < 1, is given by 

E{x, y) = \py{x, y) ■ v(x, v) ^Y.Y. CijT*{x)T*{y), (8.4) 



where 



and 



c =- C C ^^(^' y) ■ ^(^' y)T*{x)T*{y)dxdy 

Jo {Ax{l-x)y{l-y)f^ ' ^ ' ' 

/-i THxfTUy^dxdy 

Nij= / / ' T7^. (8.6) 

Jo Jo {Ax{l - x)y{l - y)f'^ 
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The denominator in these integrals is the weight function for the Chebyshev polynomials 
taking into account the change of argument given before (8.4). The first integration can 
be approximated by an SPH summation and the second can be evaluated exactly. We get 

as our estimate of the coefficients Cij. Tests on known functions shows that these estimates 

are sufficiently accurate for modes such that both i and j satisfy i < S/ {5dp) where dp is 
the intial particle spacing and S is the length of a side of the square box. 
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